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real to the complex domain — i.e., no extra complexity results when dealing 
with leaky structures or those with material/metal loss. Unlike several 
existing numerical approaches, our algorithm exhibits markedly-reduced 
sensitivity to the initial guess and allows for straightforward implementation 
on a pocket calculator 
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1. Introduction 

With an ever-improving ability to fabricate miniature functional components, the field of 
nanophotonics is in a phase of explosive growth. Not only are novel phenomena being rou- 
tinely observed in wavelength- and sub-wavelength-scale components, but these concepts are 
being rapidly translated into exotic devices for applications ranging from solar energy and in- 
formation technology, to biology IIIIIJIO. At the heart of this development is the ability of 
nanostructures to confine, guide, and scatter light in ways that are not achievable with bulk 
materials. It follows that understanding and improving the performance of these nanostructures 
requires a detailed knowledge of the electromagnetic modes they support. 

The allowed electromagnetic modes in microphotonic structures are determined by solving 
Maxwell's equations for the given geometry. The last step in this procedure is the application of 
boundary conditions at the interfaces, yielding the dispersion equations which must be solved to 
obtain the allowed modes. These dispersion equations are in general transcendental and except 
in a few simple cases, their solutions cannot be expressed in terms of elementary mathematical 
functions. 

When the solutions of the dispersion equations are known to be real, they are usually deter- 
mined using a graphical search algorithm such as the bisection method |4|. However, in lossy 
material systems or low index-contrast asymmetric waveguides, the dispersion equations have 
complex solutions [31 . This calls for a search in two dimensions (i.e., in the complex plane) and 
severely limits the effectiveness of graphical algorithms. 

Material loss and waveguide leakage are even more commonly encountered in plasmonic 
waveguides |[T]. Due to the growing technological and scientific clout of plasmonic systems 
ll6l 171 [8l l9l [TOl [TTl fT2l Q 3 1 various techniques for solving plasmon waveguide dispersion equa- 
tions are in current use. Examples include the reflection-pole method (RPM) lfT4l[T5]| . Newton's 
method fA\, and the argument-principle method fTE\. Practical implementation of these meth- 
ods requires careful programming customized for solving the problem at hand. In the methods 
that rely on curve-fitting, increasing the accuracy beyond a few decimals is often a challenging 
task, often requiring one to iterate the procedure manually based on previous results. Addi- 
tionally, the root-finding algorithms of several commercial mathematical software suites (e.g., 
Mathematica'"'^) can be very sensitive to initial guesses provided by the user; estimation of this 
initial guess is especially difficult when the solution is complex. 

In this paper, we present an easy-to-implement iterative procedure for solving complex tran- 
scendental dispersion equations that is relatively insensitive to initial guess. Our method applies 
to rectangular multilayer dielectric and plasmonic waveguides which may have either material 
or leakage loss. We first use a simple numerical example to illustrate the use of this technique 



Fig. 1. Graphs of the left- (red) and the right- (blue) hand sides of Eq. (|4]l. The red dot 
indicates the approximate location of the solution near x ~ 0.5. 



and its convergence behavior. We then successively apply the procedure to find modes of di- 
electric slab waveguides, photonic wire waveguides, and plasmonic waveguides. 

Papers in the past have discussed the use of iterative techniques for solving transcendental 
equations in general |17|. Our aim in this paper is to demonstrate how this method can be 
applied to automate the design and analysis of waveguide structures of significant contemporary 
interest. 

2. Iterative technique: A simple example 

Although the general philosophy of iterative methods is well known, we will use a simple 
example to introduce our terminology, illustrate the procedure use, and build an intuitive un- 
derstanding of convergence issues. Suppose we need to solve the equation: 

F(x)=0. (1) 

where F{x) is a combination of elementary mathematical functions. We reformulate this as an 
iterative problem by converting Eq. ([TJ to: 

x^f(x). (2) 

where / is obtained by manipulating the parent function F . Next, we start with an initial guess 
x\ and obtain a sequence according to: 

Xn+\=f{Xn)- (3) 

If the sequence defined by Eq. ([3]) converges, then the limiting value is the solution to Eq. Q. 
We will illustrate the method by a simple example. Consider the transcendental equation 

sinx = 1 — X. (4) 

Fig.[T]shows the graphs of the left-hand side (LHS) and right-hand side (RHS) of Eq. (j4]), which 
intersect around x ~ 0.5. To find the root more accurately by the iterative method, we recast the 
equation in a form similar to Eq. (j2]l by choosing f(x) = 1 — sinx, setting up an iteration scheme 
following Eq. (|3]l: 

Xn+\ = 1 — sinx„. (5) 

Fig. |2|b) plots the first fifty iterations of Eq. (j5]l and indicates that the sequence converges 
to X = 0.5109734293885691. We compare the 16-decimal agreement between this value and 
the solution x = 0.5109734293885691 computed by Mathematica™'s FindRoot function. 
Fig.|2|c) shows the LHS and RHS of the iteration scheme in Eq. (j5]l and Fig.|2|a) shows how 
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Fig. 2. Solution of the example transcendental equation via the iterative method. The graphs 
(a), (d), and (g) show the left- (magenta) and the right- (gray) hand sides of Eq. j6|, and 
(|7|. The red dot indicates the position of the initial guess. Convergence/divergence behavior 
of the real [(b), (e), (h)] and imaginary [(c), (f), (i)] parts of the iterates. 



the procedure converges to the intersection point of the two curves, starting from the initial 
guess. Making an initial guess corresponds to choosing a point on the curve y = x shown by 
the red curve in Fig.|2|a). This point is indicated in Fig.|2|a) by the red dot. The vertical lines 
in the spkal represent computation of /(x) for the chosen x and the horizontal lines represent 
re-substitution of x by f{x) for the next iteration. 

Before we proceed to apply the iterative technique to practical waveguide problems, it is 
useful to gain an intuitive understanding of how the above scheme converges to the solution. 
This is important since, for any given F{x), there are usually multiple ways to choose f{x) 
which differ in their convergence behavior. For example, an equally-legitimate way of setting 
up an iteration scheme for Eq. Q would be to choose f{x) = sin^^(l — x) and obtain the 
sequence {x„} according to: 

x„+i = sin"'(l -x„). (6) 

The behavior of successive iterates is shown in Fig.|2|d-f). Regardless of how close the initial 
guess is to the actual solution, the iterative scheme diverges (spirals away) from the intersec- 
tion point. Although the real part of the solution appears to converge after about 25 iterations, 
the overall complex solution does not converge, as seen from the undamped oscillations in the 
imaginary part of the solution. Thus, this iteration scheme, although derived from the same par- 
ent equation, does not converge. This behavior is typical for transcendental equations involving 
trigonometric functions and is therefore encountered for waveguide problems, as we will show 
in section 152] 

On the other hand, convergence of Eq. Q can be improved considerably by choosing an 
alternative iterative form as follows. Squaring both sides of Eq. Q yields; 



sin^x = 1 -cos^x = (1 — x)^ = 1 — 2x + x^. 



(7) 




Fig. 3. Criterion for convergence of tlie iterative solution. The magenta line in both figures 
is the curve y = x and the blue lines are two different cases of j = /(x). The solution (a) 
diverges for |/'(x)| > 1 and (d) converges for |/'(x)| < 1. (b), (c), (e), and (f) show how the 
convergence/divergence is reflected in the behavior of the real and the imaginary parts of 
the successive iterates. 



from which we can obtain yet another sequence of {;<;„} of iterations according to: 

Xn+l = ^ {xl + COS^ Jt:„) . (8) 

Fig. 2 g-i) depict the convergence of succesive iterates of Eq. (j8]l. By comparing Fig.|2|b) with 
Fig. 2Ti) it is apparent that Eq. (j8]l converges to the solution significantly faster than Eq. (j5]l. 

The preceding examples illustrate the crucial importance of choosing appropriate iterative 
forms. Arriving at such a form necessitates an understanding of the convergence of the itera- 
tive scheme. Figure |3] pictorially shows the criterion for convergence of the iterative scheme 
in Eq. Q. Whether the iteration spirals toward the intersection point (the solution) or away 
from it is governed by the relative slopes of the intersecting curves. If f{x) is steeper than x 
(i.e., if |/'(jc)| > 1), the successive computations and re-substitutions spiral away from the in- 
tersection and the solution diverges. On the other hand, if f{x) rises gently compared to x (i.e., 
if |/'(x)| < 1), then the scheme spirals toward the intersection and the solution converges. For 
convergent functions, the rate of convergence is governed by the magnitude of \ f'{x) \ . The RHS 
of Eq. ^ is much flatter {\f'{x) | ^ 1) near the solution compared to the RHS of Eq. (jsj), which 
results in the latter's considerably slower convergence. While the preceding justification is not a 
rigorous proof, it provides a basis for choosing the manipulations required to obtain convergent 
forms of the practically-useful equations we consider in the forthcoming sections IJ8J . 



3. Dispersion equation of a general asymmetric three-layer structure 

The design of many integrated photonics systems of practical importance — including both 
novel plasmonic waveguide architectures as well as conventional dielectric waveguides — can 
often be reduced to solving for the effective mode index and electromagnetic field distributions 
in a three-layer planar structure. As such, we focus our analysis on a generalized three-layer slab 
waveguide, as shown in Fig.|4|a). The central layer is called the core and has a complex relative 
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Fig. 4. Geometries and modes of three-layer infinite slab waveguide structures. Parts (e-i) 
plot the typical magnetic field profiles. 



permittivity £f. The bottom and top cladding layers are called the substrate (with permittivity 
Es) and cover (with permittivity Eq), respectively. 

Depending on which materials comprise the three layers, it is possible to classify the most 
commonly-encountered waveguides into three basic types. The dielectric waveguides, shown in 
Fig.|4|b), have all three layers made of dielectrics, with gf > fisjEc. The structures in Figs.Wc) 
and (d) are two basic configurations of plasmonic waveguides. The first, shown in Fig. ^c), 
consists of an dielectiic sandwiched between two (possibly different) metals and is commonly 
known as a metal-dielectric -metal (MDM) waveguide. The second structure, shown in Fig.|4|d), 
consists of a metal layer between two (possibly different) dielectrics; it is commonly known as 
an dielectric -metal-dielectric (DMD) waveguide. 

Figures |4|e-i) schematically show the possible modes of a general three-layer slab wave- 
guide. These modes are distinguished according to several criteiia: mode number (number of 
zero-crossings in the field), polarization (transverse electric or transverse magnetic), field sym- 
metry (even or odd), and field profile in the core (sinusoidal or hyperbolic). Of the possible field 
profiles shown in Fig.|4] the mode shown in (e), having sinusoidal core-field variation with no 
zero-crossings, is exclusive to dielectric waveguides and is known as the fundamental mode. 
Modes (h) and (i), having hyperbolic core-field variation, are exclusive to plasmonic waveg- 
uides and are known as the gap-plasmon modes. Modes (f) and (g), having sinusoidal core-field 
variation with > 1 zero-crossings, are common to both dielectric and plasmonic waveguides. 

This rich variety of modes aiises out of the solutions of the dispersion equation written for 
an asymmetric three-layer slab waveguide. In principle, a single dispersion equation can com- 
pletely describe all the modes of both dielectric and plasmonic waveguides. However, because 
the nominally-assumed core-field distributions for the two types are different (sinusoidal for di- 
electric waveguides and hyperbolic for plasmonic waveguides), we choose to write two separate 
equations for the two types l|T9ll20l : 

tan(kh) = — ^ ... for dielectric waveguides. (9) 



tanh( Kh) = ^ ... for plasmonic waveguides. (10) 

For convenience, we have collected the symbol definitions and their expressions in Table [T] 
Using these definitions, Eqs. (|9| and ( [TO| i can be cast explicitly as transcendental equations in 
a single complex variable k or K. Furthermore, it is evident that Eq. (|9]l transforms to Eq. ( [T0| 
for purely imaginary values of k (i.e., k iK, K e M). Having solved for k or K, the mode's 
complex effective index Weff may be calculated using relations in Table [T] 
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1 (forTE), Ef/gc (forTM) 
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1 (for TE), ef/£s (for TM) 




^0 V £f ^ ko^Sf — Eg 


Qc, Qs 


h)\/£c - £f, feoV^s - £f 


7c. 7s 






^Ki + K\ ^Ki + K^ 


^c, 




Gc, Gs 


^Jk^ + p^y}, ^Jk^ + q^y} 


S 


{pac + qas)/2 


A,B 




"eff 


\J £■{— {k/k(j)^ (sinusoidal core-fields), 
^ £f+ (K/ko)^ (hyperbolic core-fields) 



Table 1. Definitions of various quantities and their expressions in terms of the material 
parameters and the perpendicular core wavevector k or K. 
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Fig. 5. Categories of three-layer slab waveguides for the purposes of iterative solution. 



We will now begin the description of the actual iterative solution procedure. We saw in 
section [T] that the convergence of this technique depends crucially on the iteration function. 
As a result, the iterative form for determining modes with sinusoidal core-fields [Fig. Qe-g)] 
is different from the one needed for determining modes with hyperbolic core-fields [Fig.Qh- 
i)]. Therefore, we split our description into two broad categories: dielectric waveguides and 
plasmonic waveguides. These two categories are further divided depending upon the degree 
of confinement (strong or weak) for dielectric waveguides, and core/cladding type (MDM or 
DMD) for plasmonic waveguides. Figure [5] depicts the division of the waveguide modes that 
we have made for the purposes of describing our iterative solution process. In the forthcoming 
sections, we will obtain and test rapidly-convergent iterative forms for calculating the mode 
indices of these four sub-categories of three-layer asymmetric slab waveguides. These useful 
forms will be boxed for the convenience of the reader. 

4. Modes of dielectric waveguides 

4.1. Strong-confinement dielectric waveguides 

Using the half-angle identity tan(z/2) = — cotz± vT+coPz, Eq. ^ can be transformed to: 

t nu^\ {pq%Y.~k^)±G,Gs 

tm{kh2)^^ — ^ . (11) 
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Fig. 6. Iterative solution of a strong-confinement waveguide using Eq. {12}. (a) The con 



vergence plot depicting the normalized LHS (magenta) and RHS (gray) of Eq. l|12 1; f{k) 
on the v-axis of of (a) refers to the function on the right hand side. The convergence of the 
real and imaginary parts of the effective index are shown in (b) and (c), respectively. 



Inverting the tangent function in Eq. ( [TT] ) gives us the convergent iterative form for strong- 
confinement waveguides as: 



. Mn + tan 



{pqYcnYsn - kl) ± GcnG^, 
k„{pYcn+qYsn) 



Here M is an integer and the quantities Yen, Ysn, Gen, and Gs„ are related to the n'' 
through the relationships in Table [T] Although the form of Eq. 



(12) 



iterate k„ 



12 1 may seem daunting, we 



note that it is a general equation which can calculate odd and even modes of both transverse 
electric (TE) and transverse magnetic (TM) polarizations for asymmetric slab waveguides. All 
of the functions involved in Eq. (12 1 are found on scientific calculators and in mathematical 



software tools. The iterative procedure obviates any need for a graphical search and has an 
inherent ability to give arbitrary-precision solutions — a feat difficult to achieve using graphical 
interpolation. 

Before embarking on actual computation, we need to specify the parity (+ or — sign), order 
(M), and the polarization {p and q) of the desired mode. For TE polarization (p = q = 1) 
and even parity (+ sign), the effective indices of successive even TE modes may be simply 
calculated by setting M = 0,1, . . .m and iterating according to Eq. ( 12 1; this yields the TEq, 
TE2, . . . TE2m modes. For TE polarization with odd parity (— sign), setting M = 1,2, ... ,m 
yields the TEi, TE3, . . . TE2m-i modes. Note that for even modes, the mode order begins with 
M — Q, whereas for odd modes, it begins with M = I. The procedure for obtaining TM mode 
indices is identical, with proper input of {p,q) as shown in Table [T] 

To illustrate the use of the iterative method, we choose a l-/im-thick silicon-on-insulator 
(SOI) slab waveguide operating at 1550 nm as a test structure. The refractive indices of the 
silicon film, oxide substrate, and air cover are assumed to be ^/Sf = 3.50, -/e^ = 1.45, and 
— 1.00 respectively. Using the TEo mode as an example. Fig. |6] depicts how Eq. (12i 
converges to the solution. Fig. |6|b) and (c) show the real and imaginary values of the first 
twenty iterates of Eq. ( 12 1. Although the iteration is canied out in terms of k, we have chosen to 



depict the convergence in terms of the effective index Weff — P/ko since this is a more familiar 
quantity in waveguide analysis. Fig. |6|a) shows the LHS and RHS of Eq. ( [T2] l and how the 
iterative scheme converges to the solution. Notice that the slope condition mentioned in section 
[Tjis satisfied for this particular case. 

To examine the iterative technique in contrast with standard numerical techniques, we solved 
for the modes of the same structure using Mathematical^. Tablel2]compares the two solutions. 
Mathematica'^ calculations were performed by feeding Eq. ([Tl) to the FindRoot function 
with an initial guess dependent on the mode order, parity, and polarization. To obtain the initial 
guesses, we plotted the left- and right-hand sides of Eq. (Ill and made a manual estimate based 



Mode 




Parity 


M 


Iterative scheme 


MathematicaTM 


TEo 


(1,1) 


+ 





3.4347458991523551 


3.4347458991523551 


TEi 


(1,1) 




1 


3.2327892969869200 


3.2327892969869201 


TE2 


(1,1) 


+ 


1 


2.872310278807719 


2.872310278807719 


TE3 


(1,1) 




2 


2.302024617480549 


2.302024617480549 


TMo 


(ef/ec,ef/es) 


+ 





3.4165068626393461 


3.4165068626393461 


TMi 


(ef/ec,ef/es) 




1 


3.1541909024008027 


3.1541909024008027 


TM2 


(ef/ec,ef/es) 


+ 


1 


2.668932488161409 


2.668932488161409 


TM3 


(ef/ec,ef/es) 




2 


1.865243634178012 


1.865243634178012 



Table 2. Comparison of mode indices for a silicon-on-insulator (SOI) slab waveguide com- 
puted via the iterative method and the FindRoot function in Mathematica™ 



on the intersection of the graphs. Given the multiple solutions inherent in Eq. (Ill, these initial 
guesses needed to be fairly close to the actual solution for Mathematica'^ to return the correct 
solution. On the other hand, all iterative solutions in Table[T]were initiated with the same initial 
guess ki — Ks- Moreover, the convergence of the iterative technique is completely insensitive 
to the exact value of this initial guess: we obtained the same effective indices for initial guesses 
ranging from to lO^K^ (including complex values). While a precise mathematical character- 
ization of the convergence behavior is outside the scope of our paper, this exercise does suggest 
the robustness of the iterative technique with respect to the accuracy of the initial guess. 

4.1.1. Symmetric strong-confinement waveguides 

If the relative permittivities of the cover and the substrate are equal, then the structure is called 
a symmetric waveguide. Since several important photonic device structures employ symmetric 
waveguides, we now give explicit iterative forms for calculating their mode indices. These 
iterative forms arise out of the simplifications in Eq. ( [T2] ) due to the equality of Cs and Ec'. 



2 



MTT + tan"' (^p^Kl/kl-\^ 



Mn ~cot-^ p\ K}/kl-\ 



. for even mode s . (13a) 
. for odd modes. (13b) 



We also note that, like Eq. ( [T2] i, M assumes values starting from for even modes and 1 for 
odd modes. 

4.2. Weak-confinement dielectric waveguides 



Eq. (12 1 converges for a very wide range of refractive indices, wavelengths, and waveguide 
thicknesses which are commonly used in practical integrated photonics designs. It performs 
poorly, however, when applied to the weak-confinement single-mode waveguide structures rou- 
tinely fabricated out of III-V (e.g., GaAs/AlGaAs) or organic materials. In the context of mul- 
timode waveguides, the strong-confinement formula encounters convergence problems when 
used to calculate indices of modes near the waveguide cutoff. We highlight this limitation and 
its remedy through the following example. 

Consider a l-/im-thick air-clad GaAs waveguide on an Alo.1Gao.9As substrate operating at 
1550 nm. The refractive indices of GaAs and Alo.1Gao.9As are assumed to be JYf = 3.300 and 



Es — 3.256 respectively. We attempt to find the fundamental TE mode iteratively by inputting 




121. (a) The 



Fig. 7. (a-c) Iterative solution of a weak-confinement waveguide using Eq 
LHS and RHS of Eq. ([T2|; f{k) on the y-axis denotes the RHS of Eq. l[t2|. Notice the 
apparent convergence of the real part (b) and the oscillatory divergence in the imaginary 
part (c) of the effective index, (d-f) Solution of the same weak-confinement waveguide 
problem using Eq. 1 15 i. (d) The LHS and RHS of Eq. 1 15 i; f{k) on the y-axis denotes the 



RHS of Eq. \i5) . (e) and (f) show the convergence of the real and imaginary parts of the 
effective index. 



the corresponding parameters to Eq. ( 12 1 (M = Q, p = q = I, and + sign). Fig. |7];b) and (c) 



show the behavior of the real and imaginary parts of the calculated effective index for first 50 
iterates. The real part of the effective index initially diverges to reach a maximum at iteration 
number 30 and begins to converge thereafter Starting from the thirtieth iteration, the imaginary 
part oscillates between ±0.01. These oscillations do not damp out even if we increase the 
number of iterations to 10"*. Mathematica'"^, however, indicates that this structure supports TEq 
and TMo modes with indices of 3.266 and 3.263, respectively. We conclude that the iterative 
scheme in Eq. 



12 1 fails to converge for this example. 



This failure to converge can be understood by examining the behavior of the LHS and RHS 
of Eq. ( 12 1 near the intersection point, as shown in Fig.|7ja). Near the solution, \ f[{K) | > 1. As 
a result, the iterates begin to diverge away from the intersection point. In fact, l/J (fc) | continues 
to increase away from the intersection, causing a rapidly-increasing divergence. However, at 
K ~ 0.537A;o, abruptly becomes < 1, causing the real part of the iteration to converge. 

This convergence of the real part is misleading, however since the imaginary part shows un- 
damped oscillations. The behavior of the strong confinement formula in this case is similar to 
the example problem in section [T| illustrated in Fig.[2|d-f). 

In summary, the basic cause of the failure of convergence is the breakdown of the slope 
condition mentioned in section[T] This is remedied by recasting Eq. (|9]l into a form which satis- 



fies the convergence condition. To this end, we use the identity tan^ z ■ 
Eq. (|9]l as: 

A;^ ± GcGsCos(fe/z) = pqYcY?.- 



1/cos 



1 , rewriting 
(14) 



The convergent form suitable for weak confinement is obtained by solving for k contained in 



the 7c 7s term on the right hand side of Eq. ( 14 1: 



' {pqK,K,f - (Gc„Gs„)2cos2(fe„/z) + {p^q^ - l)fe4 



(15) 



pW {Ki+Ki)T2GcnGs„cos{k„h) 
where Gc„ and Gs„ are related to k„ through the relationships in Table [T] We emphasize that 
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3.26599646645606654 


3.26599646645606654 


TMo 


(ef/ec,ef/es) 


3.26338400537407312 


3.26338400537407312 



Table 3. Comparison of mode indices for an Alo.iGao gAs/GaAs/air slab waveguide com- 
puted via the iterative method and the FindRoot function in Mathematical"^ 



Eq. ( 15 I is a general equation capable of solving for near-cutoff TE and TM modes of odd and 
even parity in three-layer asymmetric slab waveguides. The choice of + or — sign in Eq. ( [T5] l 
depends on whether the cutoff occurs at the odd (+) or the even (— ) mode. When analyzing a 
given structure, it is difficult to determine a priori the parity of the cutoff mode. The prescription 
is therefore to use the strong-confinement formula until it fails to converge, continuing with the 
weak-confinement formula thereafter The domain of convergence of the strong- and weak- 
confinement formulas may be determined rigorously by setting the derivatives of the RHS of 
Eqs. ( [T2] l and ([TSj equal to 1 and solving for the corresponding value of k. 

Solving for the mode indices of the previously-considered GaAs/Alo iGao gAs waveguide 



using Eq. ( 15 1 yields the effective indices of the TEq and TMq modes in excellent agreement 
with the Mathematical'^ results, as seen from Table[3] The convergence behavior of the effective 
index is plotted in Fig. |7|d-f). The convergence of Eq. ( fTS) improves as the core-cladding 
index contrast decreases; similarly, the convergence of Eq. ( [T2| improves for increasing core- 
cladding index contrast. 

4.2.1. Symmetric wealc-confinement dielectric waveguides 

Similarly to strong-confinement waveguides, having identical cover and substrate indices con- 
siderably simplifies the iterative form for weakly-confined slab waveguides, which can be writ- 
ten as: 

k„+i = — — ... for even modes. (16a) 

Vl+/7-2tan2 {kJt/2) 

kn+i = — = ... for odd modes. (16b) 

v/l+p-2cot2(fc„V2) 

4.3. Extension to photonic wire waveguides 

Waveguides in which light is confined in two dimensions and propagates in the third dimension 
are known as photonic wire (PW) waveguides; optical fibers, as well as ridge, rib, and chan- 
nel waveguides are examples. Because they offer stronger light confinement compared to slab 
waveguides, PW waveguides are commonly employed in on-chip optical devices where a small 
size is essential for achieving several critical device specifications ( e.g., speed, low power, 
integration density, etc.). 

Although the exact field distribution in PW waveguides is best determined by numerical 
calculations, the effective index method (EIM) is remarkably successful at quickly calculating 
the mode indices for several common PW waveguide configurations 1 19 |. A knowledge of the 
effective index suffices for many important design problems. Figure |8] illustrates the two steps 
in the implementation of EIM. Figure |8ja) shows an idealized silica-clad SOI photonic wire 
waveguide having a width w and thickness h. In the first step, we disregard the confinement 
in the x-direction and solve for the effective index n' . Depending on the polarization of the 
desired mode, we will need to solve either the TE or TM equation. In the second step, we 
consider a silica-clad slab waveguide with a core index «' confined in the x-direction. Note that 



^_ Effective 
medium 

Fig. 8. (a) Typical geometry of a SOI photonic wire waveguide, (b) The two steps used in 
determining the mode index of the PW waveguide using the effective index method. The 
procedure is illustrated here for a iiy-polarized mode; steps are similar for a //,-polarized 
mode. 
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the orientation of the effective waveguide is orthogonal to that in the first step and hence the 
equation in this step is for the polarization orthogonal to that in the previous step. That is, if 
the first step uses the TE dispersion equation, then the second step uses the TM equation(and 
vice versa). The effective index n" determined in this second step is our final answer for the 
effective mode index of the original PW waveguide. 

As an example, suppose that we need to determine the effective index of the fundamental 
/i^ -polarized mode of a sub-micron PW waveguide shown in Fig. |9](a) with h ~ 450 nm and 
w = 300 nm. In the first step, we solve for the TE mode equation for the structure shown in 
Fig.|8jb). Since SOI waveguide constitutes a strong-confinement system, the first step is readily 
accomplished by iteration of Eq. ( [T2] i for the TE condition (namely, M = Q, p ^ q = 1). This 
yields the first effective index as «' — 3.073930677459340. In the second step, we consider an .it- 
confined slab with core-index n' and iterate Eq. (12 1 with the fundamental TM mode condition. 
This yields the effective index of the PW waveguide as n" = 2.652766507502340. 

To test the relative accuracy of the iterative method, we use the finite element method 
(FEM) to compute the mode index. For the waveguide with the above dimensions, FEM gives 
n" = 2.612594. It is evident that even for sub-mircron dimensions, the relative error defined 
I "iter^^FEM I /"fem ^^^^^ 1.5%. Figure |9];b) shows the variation of the relative error 
with waveguide dimensions for a fixed width-to-height aspect ratio of 1.5, from which we 
note the exponential drop in the relative error with increasing waveguide size. The reason for 
the decrease of accuracy with decreasing size lies in the increased interaction of the electro- 
magnetic fields with the corner regions of the waveguide. By decomposing a two-dimensional 




Fig. 9. (a) Photonic wire waveguide structure used for evaluation of the effective index 
method (EIM). (b) The relative error in the calculation of the mode index using EIM as a 
function of the waveguide width. Relative error is defined as | "'['tej^^pEM I /"fem- 
electric field of the fundamental mode of a SOI photonic wire waveguide calculated using 
(c) the finite element method and (d) the EIM. 



mode-solving problem into two one-dimensional problems, we effectively chose to neglect the 
interaction of the fields with the waveguide corners. This assumption becomes decreasingly 
accurate with decreasing waveguide size. Figure |9|c) and (d) show the y-component of the 
electric field, calculated using FEM and the iterative method, for the sub-micron waveguide 
considered above. Note the strong electric field at the waveguide comers in Fig.|9|c) calculated 
using FEM, and its absence in Fig. [9|d) calculated using the iterative method. In addition to 
highlighting this difference for small waveguides, this exercise illustrates the well-known fact 
that while EIM successfully calculates the effective indices, it does not guarantee satisfaction of 
the boundary conditions. It should therefore not be used to calculate fields for wavelength-sized 
structures. 

In spite of this limitation, the EIM remains an intuitive way of accomplishing designs which 
rely primarily on the knowledge of mode indices. Examples of such designs include several 
contemporary device structures such as SOI ring resonators, Bragg gratings, Mach-Zehnder 
interferometers, and directional couplers. By offering a simple means of solving the slab wave- 
guide dispersion equation, the iterative method greatly enhances the efficacy and implementa- 
tion speed of EIM. 

5. Modes of plasmonic waveguides 

The examples in the preceding sections indicate the relative simplicity of the iterative method 
when applied to dielectric waveguides. However, the true advantage of this technique emerges 
when dealing with systems having material loss (complex permittivity) or leakage loss (com- 
plex propagation constant). The modal indices are complex in both these cases and provid- 
ing graphical root-finders with a complex initial guess becomes difficult. As such, the self- 
converging behavior of the iterative algorithm becomes an invaluable asset. 

Propagation loss and waveguide leakage arise routinely in sub-wavelength-scale metallic 
waveguides. Intense research efforts are currently underway in the field of metal-based op- 
tics, also known as plasmonics (see, e.g., 1 1 1 and references contained therein). Extended 
metal structures support surface plasmon-polariton (SPP) modes that are electromagnetic waves 
strongly coupled to collective electron oscillations in the metal. To exploit the strong light local- 
ization achievable in plasmonic structures, a variety of waveguide configurations have recently 
been proposed and demonstrated [ LJ. However, because of its simplicity, the MDM geometry 
remains a canonical structure for achieving and studying strong light confinement. As such, a 
convenient technique for determining the optical modes of MDM waveguides is highly valu- 
able. 

We begin the description of the iterative solution technique for plasmonic waveguides by 
referring to the division of slab waveguide types illustrated in Fig.|4|c) and (d). We first consider 
the modes of an MDM-type plasmonic waveguide followed by a treatment of the DMD-type. 

5.7. Metal-dielectric-metal waveguides 

A MDM waveguide supports gap-plasmon and TM-like waveguide modes. Even and odd gap- 
plasmon modes are practically important, as they offer the strongest sub-wavelength field con- 
finement. However, TM-like modes are also of theoretical value and are necessary for gaining 
a complete understanding of reflection and transmission phenomena in MDM waveguides and 
antennas 1 16 |. Effective indices for both these types of modes can be conveniently determined 
using the iterative technique, as we show in the following. 



5.1.1. Gap-plasmon modes 

We start with our master plasmonic Eq. ( [T0| and rearrange it as: 

+ 2SKC0th{Kh) + pqUcUs — 0. 



(17) 



Here, ac,(Xs, and S are related to K as specified in Table[T] Treating this as a quadratic equation 
in K and solving yields the convergent iterative form for the MDM gap-plasmon modes as: 



ffn+l = — 5„coth(fc„/i) ± a/ S^coth {K„h) — pqa^^nU, 



(18) 



This simple-looking equation can calculate even (+ sign) and odd (— sign) gap-plasmon mode 
indices of a wide variety of deep sub-wavelength asymmetric plasmonic waveguides. We il- 
lustrate the use of Eq. ( 18 i through two examples. Our first structure is a 50-nm-thick gold- 



silica-silver slab waveguide operating at 1550 nm. The relative permittivities of gold, silica, 
and silver are assumed to be —95.92 — /10.97, 2.1025, and —143.49 — /9.52, respectively. The 
effective index of the even gap-plasmon mode, calculated using the + sign in Eq. ( 18 1, is given 
in Table|4] Fig.[TO]shows plots of the LHS and RHS of Eq. ( fTSj ) and the convergence of the real 
and imaginary parts of the effective index. 

The thin 50-nm slab considered above does not support an odd (antisymmetric) gap-plasmon 
mode. However, increasing the silica thickness to 3 fim allows the structure to support gap- 



plasmon modes of both symmetries. We calculate the effective indices iteratively using Eq. ( 18 i 
with + and — signs for odd and even gap-plasmon modes respectively. Once again, the iterative 
method agrees with the calculations of the FindRoot function in Mathematical'^ (Table|4|. 

5. 1 .2. TM-like waveguide modes 

The other class of modes supported by a MDM waveguide are the TM-like waveguide modes 
with profiles as shown in Fig.|4|f) and (g). For these modes, K is purely imaginary and plasmon 
Eq. ( [T0| transforms to dielectric Eq. (|9]l. The iterative form for obtaining the indices of these 
modes is identical to Eq. ([12]), with only a few slight differences as regards its implementation. 
For the case of dielectric waveguides, the mode index M assumed integer values starting from 
for even modes and 1 for odd modes. For MDM waveguides, this is reversed: M assumes 
integer values starting from 1 for even modes and for odd modes. This is a consequence of 
the signs of p and q being reversed for MDM waveguides due to the negative permittivity of 
the metal "claddings." Additionally, because of the large index difference between metals and 
dielectrics, the TM-like modes of MDM waveguides are almost always calculated using the 
strong-confinement formula in Eq. 



121, 



To show how TM-like modes are calculated for MDM waveguides, we choose a 300-nm 
thick silica layer sandwiched between semi-infinite gold and silver layers. The permittivities of 
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Fig. 10. (a) Normalized left- and right-hand sides of Eq. 1 18 1; /(k) refers to the RHS. 
Convergence of the real (b) and the imaginary (c) parts of the effective index for the funda- 
mental gap-plasmon mode. 
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Mathematica^M 
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H-fleld 



Even gap 
plasmon 

Odd gap 
plasmon 



Eq. 

+ 

Eq. 



18i 



(18i 



1.467915033129527- 
/0.001514007231254 
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/0.001440093524486 
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/ 1.98 1855964604849 
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Eq. 
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Eq. 
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([12|| 0.001924784371747- 
1.+ (4.90109582884017 



([12|| 0.000214216445512- 
1.- ;7. 58348752253 199 



0.001924784371747- 
/4.90109582884017 

0.000214216445512- 
/7.58348752253199 



Ag T Au 
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TM4 



Eq. (12 1 0.00592749529203- 



M = 2, 



/10.22010371292752 



0.00592749529203 -f- 
/ 10.220 1037 1292752 



Ag II AU 



Eq. ([12|| 0.01577537648440-F 
M = 2,- il2.83 149770403419 



0.01577537648440 
/12.83149770403419 



Table 4. Effective indices of various modes of MDM waveguides obtained using the itera- 
tive method, with a comparison to the direct solutions calculated using Mathematica'''". 



each material are the same as assumed in the previous subsection. The indices of the first five 
TM-like modes, calculated using the strong-confinement formula Eq. ( 12 1, are listed in Table|4] 
alongside the corresponding solutions from Mathematica''^. 

5.1.3. Symmetric MDM waveguides 

Many contemporary high-confinement architectures employ the symmetric MDM waveguide as 
their skeleton structure. Equality of the relative permittivities of the substrate and cover further 



simplify Eq. ( 18 1 for the gap-plasmon modes. We can express the resulting iterative forms as: 



K„+i = —py + K^tanh{K„h/2) ... for even gap plasmon. (19a) 



K„+i = —p^ + K^coth{K„h/2) ... for odd gap plasmon. (19b) 

5.2. Dielectric-metal-dielectric waveguides 

The second type of plasmonic waveguide structure commonly encountered in practice is an 
DMD waveguide, as shown in Fig. |4|d). In theory, metal films of arbitrary thickness in a ho- 
mogenous dielectric medium (including air) or in Kretschmann-type coupling configurations 
ETI are examples of DMD waveguides. In practice, we refer to metallic waveguides as DMD- 
type only if the SPP modes on the two metal-dielectric interfaces are coupled. Because of the 
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Fig. 11. Convergence of the real and imaginary parts of K, A, and B in Eq. (j2TJ for the case 
of a 50 nm thick silver-silica-silver waveguide operating at 1550 nm. 



rapid decay of the fields (with distance) inside metals, such a mode-coupling is possible only 
for thin (h < 100 nm) metal films. 

To obtain the iterative form for determining mode indices of DMD waveguides, we write 



Eq. (lOi as: 

—2Ak 

tmihKh= ^2+a2_b2 ^ (20) 

where A and B are defined in Table [T] Considering this as a quadratic equation in A and B 
leads us to the desired iterative forms. We start by identifying the lower-index dielectric as 
the substrate and using initial guesses of Aq — ko, Bq — 0, and Kq = ko- We then iterate to 
obtain five different quantities successively, using the following equations in the order shown: 



a*^-K„ coth Knh±J Bl + K^cficY? (Knh), (21a) 



K = y af + kI + 2a* K„ coth^{K„h) , (21b) 
K-„+i = ^{^^^+b*fJp^+Ql (21c) 

An+\ = {P^c.n+l + q^i.n+l) /2, (21d) 
Bn+l = O'^c.n+l -?^s,n+l)/2. (21e) 

Here, a*,b* are intermediate dummy parameters and £,c.„+i and £,s,„+i are related to (Cn+i 
through Table[T] The indices of the low- and high-energy plasmon modes are obtained by using 
the + and the — signs, respectively, in Eq. ( |21a[ ). The large asymmetry in the decay constants 
in the metal core (k) and dielectric claddings (^c,^s) makes the determination of the effective 
indices of DMD modes a numerically-challenging problem. For the iterative procedure, this 
translates into a difficulty in obtaining convergent iteration expressions. Therefore, unlike the 
previous cases, the iterative procedure for DMD waveguides involves multiple steps instead of 
a single closed-form iteration function. Admittedly, this multi-step algorithm goes against the 
"pocket calculator" philosophy, but nevertheless is convenient and efficient compared to direct 
numerical solutions. 

We will now illustrate the use of Eq. ( [2T] l through specific examples. Consider the operation 
of a 50-nm-thick silver film on a silica substrate at 1550 nm. The relative permittivities of 
silver and silica are the same as assumed previously; the top cover is air (e^ = 1). We start 
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Table 5. Effective indices of various modes of an DMD waveguide obtained using the 
iterative method and their comparison with direct solutions using Mathematical"^. 



by specifying Aq^Bq, and fCo, and obtain their successive values according to Eq. (21 1. The 
convergence of A„ and K„ is shown in Fig. [TT[ Table |5] compares the solution obtained 
using the iterative scheme with the direct solution computed by the Mathematical^ FindRoot 
function. In our computations, the FindRoot function was unable to return the complex mode 
index unless we gave a very precise initial guess for both the real and the imaginary parts. On 
the other hand, the iterative method computed the complex mode index regardless of the initial 
guess. 

Our next example is a symmetric structure with a 100-nm silver film sandwiched between 
two silica layers. Although 100 nm is at the boundary of the film behaving like an DMD wave- 
guide versus two separate interfaces, we choose these dimensions to highlight the robustness 
of the iterative technique even for the most extreme cases of root-finding. For this structure, 
both high- and low-energy modes exist, whose indices are conveniently calculated by using the 



— and + signs, respectively, in Eq. (21 1. The indices calculated using the iterative method and 
Mathematical^ foj- fjjjg example are included in Table [5] and show good agreement. Calculat- 
ing these indices was especially difficult using FindRoot because of their close numerical 
proximity. The iterative technique, on the other hand, specifies separate functions (the + and— 
signs) which are guaranteed to converge to these different modes, irrespective of their numerical 
proximity. 

5.2.1. Symmetric DMD waveguides 

Symmetric DMD waveguides appear in the form of idealized waveguide geometries such as 
metal films in homogenous dielectric media (including air). Fortunately, for the symmetric 
case, the iterative scheme consists of a single equation which can be written as: 

2c 

(Cn+i = — . . . for even plasmon mode. (22a) 

y^l-/7-2tanh2(K-„/!/2) 

Qc 

K„+i = — ... for odd plasmon mode. (22b) 



'l-/7-2coth2(K-„/i/2) 

These equations do not follow automatically from Eq. ( [2T| but instead have to be derived sepa- 



rately by considering the dispersion equation for the symmetric DMD waveguide. Eq. (22 1 are 



preferable to Eq. (21 1 for analyzing symmetric DMD structures owing to better convergence 



behavior and an obvious ease in programming. 



6. Conclusion 



We have presented a robust and an easy-to-implement iterative technique for determining com- 
plex propagation constants of asymmetric dielectric and plasmonic waveguides. At the heart of 
our procedure are the closed-form iteration functions, namely the boxed Eqs. ([T2|, ([Tsj, ( [Ts] ), 
and ( |2T] i. In addition to the programming ease, the iterative technique has an inherent ability 
to give arbitrary-precision answers — a feat difficult to achieve using graphical or curve-fitting 
algorithms such as the reflection-pole method. Because of its insensitivity to initial guess, it is 
our hope that this technique will help facilitate design tasks that require rapid and automated 
calculation of mode indices for a variety of nanophotonic structures having a rectangular ge- 
ometry. 

Acknowledgements 

We thank Dr. E. Kocaba§ and Prof. D. A. B. Miller for several helpful discussions. We also 
gratefully acknowledge the financial support from the Si-based Laser Initiative of the Multi- 
disciplinary University Research Initiative (MURI) under the Air Force Aerospace Research 
Award No. FA9550-06- 1-0470 and the MARCO Interconnect Focus Center. 



